Generalized Procrustes Analysis

• Align objects to common coordinate system using superimposition (GPA)

\[\small\mathbf{Z}=\frac{1}{CS}\mathbf{(Y-\overline{Y})H}\]

• Standardize ‘nuisance’ parameters (\(\small{CS}\), \(\small\mathbf{H}\), \(\small\mathbf{\overline{Y}}\)) as:
  • \(\small\mathbf{\overline{Y}}\) = translate centroid to the origin
  • \(\small{CS}\) = scale by centroid size
  • \(\small\mathbf{H}\) rigid rotation (based on optimal angle using all landmarks)

• Works for both 2-D and 3-D landmarks data
• THE fundamental approach for morphometrics

The Problem with Curves I

• Landmarks provide rigorous quantification of shape
• Some objects have few landmarks, but have defined boundaries (curves) or surfaces

• How do we capture shape information from such structures?

The Problem with Curves II

• Some objects contain multiple ‘types’ of morphometric data
• Types typically analyzed separately

• How do we quantify the overall shape of such structures?

Mathematically-Defined Points: Problems

• Define landmark homology mathematically (e.g., equal spacing)
• Often produce suggestive (but incorrect) results

Mathematically-Defined Points: Problems

• Define landmark homology mathematically (e.g., equal spacing)
• Often produce suggestive (but incorrect) results

• Implied shape differences:

Mathematically-Defined Points: Problems

• Define landmark homology mathematically (e.g., equal spacing)
• Often produce suggestive (but incorrect) results

• Implied shape differences:

• Can we estimate the shape difference without adding the twisting?

Example from: Gunz et al. (2005) in Modern morphometrics in physical anthropology.

Mathematically-Defined Points: Problems (Cont.)

• Another example

Mathematically-Defined Points: Problems (Cont.)

• Another example

• Implied shape differences:

Mathematically-Defined Points: Problems (Cont.)

• Another example

• Implied shape differences:

• Can we estimate the shape difference without adding the twisting?

Example from: Gunz et al. (2005) in Modern morphometrics in physical anthropology.

Semilandmarks and Procrustes Relaxation

• Procrustes superimposition of 2D points and curves (Bookstein 1997)
• Points along curves are ‘degenerate’ landmarks (semilandmarks)
• Semilandmarks constrained to be on the curve (thus have \(\small{\approx{1}}\) df)

Semilandmarks and Procrustes Relaxation

• Procrustes superimposition of 2D points and curves (Bookstein 1997)
• Points along curves are ‘degenerate’ landmarks (semilandmarks)
• Semilandmarks constrained to be on the curve (thus have \(\small{\approx{1}}\) df)
• Digitize a set number of points along an outline, and allow them to ‘slide’ along curve to improve the fit of a specimen to a reference curve: minimize difference between them

Semilandmarks and Procrustes Relaxation

• Procrustes superimposition of 2D points and curves (Bookstein 1997)
• Points along curves are ‘degenerate’ landmarks (semilandmarks)
• Semilandmarks constrained to be on the curve (thus have \(\small{\approx{1}}\) df)
• Digitize a set number of points along an outline, and allow them to ‘slide’ along curve to improve the fit of a specimen to a reference curve: minimize difference between them
• Treat as points for standard landmark analysis (GPA+TPS)

Bookstein (1997) Med. Image Anal.

Procrustes Relaxation: Conceptual Procedure

• Incorporate ‘sliding’ of semilandmarks into GPA algorithm

• GPA algorithm with sliding landmarks: Procedure

  • 1: Translate and scale specimens

Procrustes Relaxation: Conceptual Procedure

• Incorporate ‘sliding’ of semilandmarks into GPA algorithm

• GPA algorithm with sliding landmarks: Procedure

  • 1: Translate and scale specimens
  • 2: LS rotation, find consensus

Procrustes Relaxation: Conceptual Procedure

• Incorporate ‘sliding’ of semilandmarks into GPA algorithm

• GPA algorithm with sliding landmarks: Procedure

  • 1: Translate and scale specimens
  • 2: LS rotation, find consensus
  • 3: ‘Slide’ landmarks along curve in some way

Procrustes Relaxation: Conceptual Procedure

• Incorporate ‘sliding’ of semilandmarks into GPA algorithm

• GPA algorithm with sliding landmarks: Procedure

  • 1: Translate and scale specimens
  • 2: LS rotation, find consensus
  • 3: ‘Slide’ landmarks along curve in some way
  • 4: Perform GPA to obtain new consensus

Procrustes Relaxation: Conceptual Procedure

• Incorporate ‘sliding’ of semilandmarks into GPA algorithm

• GPA algorithm with sliding landmarks: Procedure

  • 1: Translate and scale specimens
  • 2: LS rotation, find consensus
  • 3: ‘Slide’ landmarks along curve in some way
  • 4: Perform GPA to obtain new consensus
  • 5: Iterate steps 2-5 to improve fit

Procrustes Relaxation: Conceptual Procedure

• Incorporate ‘sliding’ of semilandmarks into GPA algorithm

• GPA algorithm with sliding landmarks: Procedure

  • 1: Translate and scale specimens
  • 2: LS rotation, find consensus
  • 3: ‘Slide’ landmarks along curve in some way
  • 4: Perform GPA to obtain new consensus
  • 5: Iterate steps 2-5 to improve fit

• Two analytical issues:

  • What direction to slide?
  • How far to slide?

What Direction to Slide?

• Slide landmarks along curve

• Adjacent points on curve define sliding directions: \(\small{U}_{ij}=[(x_{j-1}-x_{j+1}),(y_{j-1}-y_{j+1})]\)

What Direction to Slide?

• Slide landmarks along curve

• Adjacent points on curve define sliding directions: \(\small{U}_{ij}=[(x_{j-1}-x_{j+1}),(y_{j-1}-y_{j+1})]\)

• Assemble into tangent direction matrix

\[\tiny\mathbf{U}=\begin{bmatrix} 0 & 0 & 0 \\ u_{2x} & 0 & 0\\ 0 & u_{3x} & 0 \\ 0 & 0 & u_{4x} \\ 0 & 0 & 0 \\ \hline 0 & 0 & 0 \\ u_{2y} & 0 & 0\\ 0 & u_{3y} & 0 \\ 0 & 0 & u_{4y} \\ 0 & 0 & 0 \\ \end{bmatrix}\]

What Direction to Slide?

• Slide landmarks along curve

• Adjacent points on curve define sliding directions: \(\small{U}_{ij}=[(x_{j-1}-x_{j+1}),(y_{j-1}-y_{j+1})]\)

• Assemble into tangent direction matrix

\[\tiny\mathbf{U}=\begin{bmatrix} 0 & 0 & 0 \\ u_{2x} & 0 & 0\\ 0 & u_{3x} & 0 \\ 0 & 0 & u_{4x} \\ 0 & 0 & 0 \\ \hline 0 & 0 & 0 \\ u_{2y} & 0 & 0\\ 0 & u_{3y} & 0 \\ 0 & 0 & u_{4y} \\ 0 & 0 & 0 \end{bmatrix}\]

• Landmark positions after sliding found as: \(\small\mathbf{Y}^*=\mathbf{Y^0-UT}\)

How Far to Slide?

• Most conservative approach: slide to minimize shape differences
• Sliding defined by \(\small\mathbf{T}\): \(\small\mathbf{Y}^*=\mathbf{Y^0-UT}\)
• Two criteria: bending energy (\(\small{BE}\)) & Procrustes distance (\(\small{D}_{Proc}\))

How Far to Slide?

• Most conservative approach: slide to minimize shape differences
• Sliding defined by \(\small\mathbf{T}\): \(\small\mathbf{Y}^*=\mathbf{Y^0-UT}\)
• Two criteria: bending energy (\(\small{BE}\)) & Procrustes distance (\(\small{D}_{Proc}\))
• 1: Slide landmarks to minimize \(\small{BE}\) via GLS: \[\small\mathbf{T}=(\mathbf{U^T\mathcal{L}_p^{-1}U})^{-1}\mathbf{\mathcal{L}_p^{-1}U^T}(\mathbf{Y^0-Y_{ref}})\]

where \(\small\mathbf{\mathcal{L}_p^{-1}}=\begin{bmatrix} \mathbf{L}_p^{-1} & 0 \\ 0 & \mathbf{L}_p^{-1} \end{bmatrix}\)

How Far to Slide?

• Most conservative approach: slide to minimize shape differences
• Sliding defined by \(\small\mathbf{T}\): \(\small\mathbf{Y}^*=\mathbf{Y^0-UT}\)
• Two criteria: bending energy (\(\small{BE}\)) & Procrustes distance (\(\small{D}_{Proc}\))
• 1: Slide landmarks to minimize \(\small{BE}\) via GLS: \[\small\mathbf{T}=(\mathbf{U^T\mathcal{L}_p^{-1}U})^{-1}\mathbf{\mathcal{L}_p^{-1}U^T}(\mathbf{Y^0-Y_{ref}})\]

where \(\small\mathbf{\mathcal{L}_p^{-1}}=\begin{bmatrix} \mathbf{L}_p^{-1} & 0 \\ 0 & \mathbf{L}_p^{-1} \end{bmatrix}\)

• 2: Slide landmarks to minimize \(\small{D}_{Proc}\): \[\small\mathbf{T}=\mathbf{U^T}(\mathbf{Y^0-Y_{ref}})\]

Note: \(\small\mathbf{Y^0}\) and \(\small\mathbf{Y_{ref}}\) are assembled by \(\small{x}\) then \(\small{y}\) coordinates

Statistical note: Sliding by \(\small{D}_{Proc}\) assumes independence of semilandmarks, as \(\small(\mathbf{U^T\mathcal{L}_p^{-1}U})^{-1}\mathcal{L}_p^{-1}=\mathbf{I}\) (this is not realistic).

Procrustes Relaxation: Flow of Computations

• Sliding incorporated into GPA algorithm

Procrustes Relaxation: Flow of Computations

• Sliding incorporated into GPA algorithm
  • 1: Translate and scale specimens
  • 2: LS rotation, obtain reference and \(\small\mathbf{L}_p^{-1}\)

Procrustes Relaxation: Flow of Computations

• Sliding incorporated into GPA algorithm
  • 1: Translate and scale specimens
  • 2: LS rotation, obtain reference and \(\small\mathbf{L}_p^{-1}\)
  • 3: Calculate tangent directions \(\small\mathbf{U}\) \[\small{U}_{ij}=[(x_{j-1}-x_{j+1}),(y_{j-1}-y_{j+1})]\]

Procrustes Relaxation: Flow of Computations

• Sliding incorporated into GPA algorithm
  • 1: Translate and scale specimens
  • 2: LS rotation, obtain reference and \(\small\mathbf{L}_p^{-1}\)
  • 3: Calculate tangent directions \(\small\mathbf{U}\) \[\small{U}_{ij}=[(x_{j-1}-x_{j+1}),(y_{j-1}-y_{j+1})]\]
  • 4: Slide landmarks along tangents using \(\small{BE}\) or \(\small{D}_{Proc}\) \[\small\mathbf{T}=(\mathbf{U^T\mathcal{L}_p^{-1}U})^{-1}\mathbf{\mathcal{L}_p^{-1}U^T}(\mathbf{Y^0-Y_{ref}})\] \[\small\mathbf{T}=\mathbf{U^T}(\mathbf{Y^0-Y_{ref}})\]

Procrustes Relaxation: Flow of Computations

• Sliding incorporated into GPA algorithm
  • 1: Translate and scale specimens
  • 2: LS rotation, obtain reference and \(\small\mathbf{L}_p^{-1}\)
  • 3: Calculate tangent directions \(\small\mathbf{U}\) \[\small{U}_{ij}=[(x_{j-1}-x_{j+1}),(y_{j-1}-y_{j+1})]\]
  • 4: Slide landmarks along tangents using \(\small{BE}\) or \(\small{D}_{Proc}\) \[\small\mathbf{T}=(\mathbf{U^T\mathcal{L}_p^{-1}U})^{-1}\mathbf{\mathcal{L}_p^{-1}U^T}(\mathbf{Y^0-Y_{ref}})\] \[\small\mathbf{T}=\mathbf{U^T}(\mathbf{Y^0-Y_{ref}})\]
  • 5: Perform GPA and obtain new reference
  • 6: Iterate 2 - 5 to improve fit

Semilandmarks: Differing Approaches

• Visual differences between \(\small{BE}\) and \(\small{D}_{Proc}\) sliding:

Mathematical comments: Using \(\small{BE}\) is akin to a weighted regression, which accounts for spatial proximity (non-independence) of adjacent landmarks. This is a more reasonable model, as one uses semilandmarks precisely to represent curves by an adjacent set of points which are explicitly not independent

Procrustes Relaxation: Refinements

• Semilandmarks slide along tangent directions
• But if slid to far they ‘leave’ the actual curve
• Solution is to project back to curve during GPA + sliding

Combining Points and Curves

• One can allow landmarks to remain ‘fixed’ while semilandmarks slide (e.g. endpoints above)
• Information incorporated into \(\small\mathbf{U}\)
• Columns of \(\small{0}\) represent landmarks that don’t slide

Combining Points and Curves

• One can allow landmarks to remain ‘fixed’ while semilandmarks slide (e.g., endpoints above)
• Information incorporated into \(\small\mathbf{U}\)
• Columns of \(\small{0}\) represent landmarks that don’t slide

\[\tiny\mathbf{U}=\begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & u_{2x} & 0 & 0& 0 \\ 0 & 0 & u_{3x} & 0 & 0 \\ 0 & 0 & 0 & u_{4x} & 0 \\ 0 & 0 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 & 0 & 0 \\ 0 & u_{2y} & 0 & 0 & 0\\ 0 & 0 & u_{3y} & 0 & 0 \\ 0 & 0 & 0 & u_{4y} & 0 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix}\]

Note: with this notation, generalized inverses are required to complete GLS sliding (and process can be SLOW)

Note 2: geomorph uses very clever algebra to eliminate redundancies in computations to speed this process up considerably

Combining Landmarks and Curves: Example

• Cranial profiles in Homo

• External differences in frontal bones of archaic and modern humans are well-known

• Archaic humans have ‘frontal flattening’ while modern human crania are more vertically rounded

• Common anthropological interpretation: more-rounded crania in modern humans due to enlargement of frontal lobes of brain

• Profiles of internal and external frontal bones were digitized from 5 mid-Pleistocene and Neanderthal crania, and 16 modern humans and compared

Bookstein et al. (1999). Anatomical Record.

Semilandmark Example: The Data

Semilandmark Example: The Data

Semilandmark Example: Results

• Large differences in external profiles, but NO difference of internal profiles

Semilandmarks on 3D Curves

• Semilandmarks on 3D curves work identically, by defining tangent directions of \(\small\Delta{x}\), \(\small\Delta{y}\), and \(\small\Delta{z}\)

\[\tiny\mathbf{U}=\begin{bmatrix} 0 & 0 & 0 \\ u_{2x} & 0 & 0\\ 0 & u_{3x} & 0 \\ 0 & 0 & u_{4x} \\ 0 & 0 & 0 \\ \hline 0 & 0 & 0 \\ u_{2y} & 0 & 0\\ 0 & u_{3y} & 0 \\ 0 & 0 & u_{4y} \\ 0 & 0 & 0 \\ \hline 0 & 0 & 0 \\ u_{2z} & 0 & 0\\ 0 & u_{3z} & 0 \\ 0 & 0 & u_{4z} \\ 0 & 0 & 0 \end{bmatrix}\]

• Slide as before:

\[\tiny\mathbf{T}=(\mathbf{U^T\mathcal{L}_p^{-1}U})^{-1}\mathbf{\mathcal{L}_p^{-1}U^T}(\mathbf{Y^0-Y_{ref}})\] \[\tiny\mathbf{\mathcal{L}_p^{-1}}=\begin{bmatrix} \mathbf{L}_p^{-1} & 0 & 0 \\ 0 & \mathbf{L}_p^{-1} & 0 \\ 0 & 0 & \mathbf{L}_p^{-1} \end{bmatrix}\]

\[\tiny\mathbf{T}=\mathbf{U^T}(\mathbf{Y^0-Y_{ref}})\]

Semilandmarks on Surfaces

• 3D surface landmarks slide in the plane
• Plane defined by PC1 & PC2 of surrounding points

• \(\small\mathbf{U}\) has 2 sets of tangent directions Slide using \(\small{BE}\) or \(\small{D}_{Proc}\) as before

\[\tiny\mathbf{U}= \left[ \begin{array}{ccc|ccc} 0 & 0 & 0 & 0 & 0 & 0 \\ u_{2x} & 0 & 0 & w_{2x} & 0 & 0 \\ 0 & u_{3x} & 0 & 0 & w_{3x} & 0 \\ 0 & 0 & u_{4x} & 0 & 0 & w_{4x} \\ 0 & 0 & 0 & 0 & 0 & 0\\ \hline 0 & 0 & 0 & 0 & 0 & 0\\ u_{2y} & 0 & 0 & w_{2y} & 0 & 0 \\ 0 & u_{3y} & 0 & 0 & w_{3y} & 0 \\ 0 & 0 & u_{4y} & 0 & 0 & w_{4y} \\ 0 & 0 & 0 & 0 & 0 & 0\\ \hline 0 & 0 & 0 & 0 & 0 & 0\\ u_{2z} & 0 & 0 & w_{2z} & 0 & 0\\ 0 & u_{3z} & 0 & 0 & w_{3z} & 0 \\ 0 & 0 & u_{4z} & 0 & 0 & w_{4z} \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array} \right]\]

NOTE: One can see that these matrices are becoming VERY large. That will lead to long computation time (though recall geomorph uses clever algorithms to speed this up! [ask Mike!])

Defining Semilandmarks on Surfaces

• Surface semilandmarks must retain positional correspondence
• How to digitize?
• 1: ‘By hand’ (digitize same set of points from all surfaces)

Defining Semilandmarks on Surfaces

• Surface semilandmarks must retain positional correspondence, so how to digitize?
• 1: ‘By hand’ (digitize same set of points from all surfaces)
• 2: Semi-automate
  • Digitize many landmarks on surface
  • ‘Thin’ to a reasonable number (150-200) on 1st specimen
  • Match these to surface points of next specimen (using TPS)
  • Repeat

Note: resampling methods required for significance testing as \(\small{p>>n}\)

Combining Points, Curves, and Surfaces

• Landmarks and semilandmarks can all be combined
  • Digitize points
  • Digitize curve landmarks
  • Digitize surface landmarks
• Define tangent directions and generate U
• GPA + Sliding

Complications: Shape Variables and Dimensionality

• Sliding procedure results in redundant shape dimensions
  • Minimum: ~1 df lost per 2D/3D curve semilandmark; and ~2 df lost per 3D surface semilandmark (e.g., ~ \(\small{2p-4-m}\) shape variables in 2D for \(\small{m}\) = # semilandmarks)
  • However, usually many more dimensions lost, because adjacent semilandmarks encode redundant shape information
• Causes MAJOR issues for statistical hypothesis testing using Procrustes residuals

Complications: Shape Variables and Dimensionality

• Sliding procedure results in redundant shape dimensions
  • Minimum: ~1 df lost per 2D/3D curve semilandmark; and ~2 df lost per 3D surface semilandmark (e.g., ~ \(\small{2p-4-m}\) shape variables in 2D for \(\small{m}\) = # semilandmarks)
  • However, usually many more dimensions lost, because adjacent semilandmarks encode redundant shape information
• Causes MAJOR issues for statistical hypothesis testing using Procrustes residuals
• Solutions:
  • 1: PCA to ‘spin out’ redundant dimensions (remove dimensions with 0.0% variance)
  • 2: Use generalized inverses in computations
  • 3: Use resampling procedures (RRPP) for significance testing

Note: #3 also required when # shape variables exceeds # specimens

Semilandmarks: Snake Example

• Examined head shape in snakes (83 species)
• Focused on aquatic species in different families
• Is there evidence of head shape convergence?

Semilandmarks: Homo Example

• Examined ‘Zuttiyeh’ fossil morphology as compared to other early hominids

Semilandmarks: Scallop Example

• Scallop species display differing life habits
• What patterns are observed at the macroevolutionary level?

Semilandmarks: Scallop Example (Cont.)

• Recessors display unique morphology and directional evolution

Semilandmarks: Scallop Example (Cont.)

• Gliders occupy morphological convergence

The Procrustes Paradigm

• With landmarks and semilandmarks we arrive at a general solution for quantifying morphology
• Quantitative information representing:
  • Discrete anatomical points
  • Curves (outlines) of structures
  • Surfaces of structures
• Can all be combined in a single analysis of shape variation